Discriminating dynamical from additive noise 
in the Van der Pol oscillator 



C. Degli Esposti Boschi a , G.J. Ortega and E. Louis 

O ' 

a Departamento de Fisica Aplicada and Unidad Asociada of the Consejo Superior 
C^l \ de Investigaciones Cientificas, 

Facultad de Ciencias, Universidad de Alicante, Apartado 99, E-03080, Alicante, 

Spain 

h Centro de Estudios e Investigaciones, Universidad National de Quilmes, 
R. S. Pena 180, 1876, Bernal, Argentine 



Q 

U 
d 



Abstract 



We address the distinction between dynamical and additive noise in time series 

analysis by making a joint evaluation of both the statistical continuity of the series 

and the statistical differentiability of the reconstructed measure. Low levels of the 

latter and high levels of the former indicate the presence of dynamical noise only, 

while low values of the two are observed as soon as additive noise contaminates the 

signal. The method is presented through the example of the Van der Pol oscillator, 

' but is expected to be of general validity for continuous-time systems. 
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Experimental time series are always blurred with additive noise coming from a 
variety of sources [1,2]. Additive noise further complicates time series analysis, 
leading to ambiguous interpretations of the basic quantities which characterise 
the dynamics, like correlation dimension and Lyapunov exponents [3,4]. In 
particular it may lead unreachable the goal of differentiating deterministic 
from stochastic dynamics [5—11]. 

Here an attempt is made to discriminate additive noise (AN) from the dy- 
namical noise (DN) the system may have. Note that the term additive used 
here denotes a component of noise that is superimposed to the underlying 



Preprint submitted to Physica D 



8 February 2008 



"clean" signal, due to the measurement process. Common synonyms found in 
the literature are measurement or observational noise. It should not be con- 
fused with the term used in the context of stochastic processes (see, for istance, 
[12]) where it often indicates a coordinate-independent stochastic forcing as 
opposite to multiplicative noise, for which the amplitude of the random terms 
depends on the coordinates themselves. Here we do not focus on this last dis- 
tinction, and we will denote generically as dynamical every noisy mechanism 
intrinsic to the system. 

Our method is based upon a fundamental, topological, property of determin- 
istic systems, namely, the differentiability of the measure along the recon- 
structed trajectory [13,10,11] or, more specifically, the continuity of its log- 
arithmic derivative. Starting from this basis we show that noise (additive or 
dynamical) decrease the differentiability of the measure. This would in princi- 
ple hinder the possibility of differentiating the two types of noise. Then we look 
at the same property (continuity) of the coordinate. Now, while AN destroys 
continuity of the coordinate, the latter is scarcely affected by DN. Thus, a low 
continuity of both the coordinate and (the log derivative of) the reconstructed 
measure would indicate the presence of AN, while if continuity is high for the 
coordinate and low for the measure the system contains DN only (see follow- 
ing section). The method, however, doesn't differentiate sharply a system with 
the two types of noise from one having solely AN. Notwithstanding, we believe 
that the present approach is a significant step forward in the understanding 
of determinism and stochasticity. To illustrate how the method works we con- 
sider the simple case of the Van der Pol oscillator [14]. No reasons are forseen 
that may indicate that the applicability of the method is limited to simple, 
non chaotic, dynamical systems, such as the one investigated here. 

The rest of the paper is organised as follows. In section 2 we make some 
general considerations on the effects of additive and dynamical noise. In section 
3 we briefly discuss the numerical methods followed to solve the dynamical 
equations with noise and the embedding process. The procedure followed to 
evaluate the (statistical) continuity is also summarised in that section. The 
results are discussed in section 4, while section 5 is devoted to the conclusions 
of our work. 



2 General considerations 

The method of the reconstructed measure along the trajectory [10,11] has the 
capability of revealing the presence of DN in a given time series by looking 
at the degree of (statistical) continuity of the logarithmic derivative of such a 
measure. To this end, it is to be considered as complementary to other tech- 
niques, essentially based on short-time predictability [5] or on smoothness in 
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phase-space [6-8]. As far as AN is concerned, in ref. [15] it is argued that meth- 
ods of the second type can be useful even in real experimental series affected 
by measurement components. A different line is that followed by Barahona 
and Poon [16], who are able to cope with rather large amounts of AN, at the 
"price" of building the analysis upon a given class of Volterra-type nonlinear 
models adjusted on the time series itself. Here, we do not specify any a priori 
dynamics and pursue the extension of the method of Ortega and Louis [10,11] 
to deal with AN. The choice is motivated by very recent results [17] which al- 
low us to believe that this method is suitable even for high-dimensional chaotic 
systems, which somehow fool the above-mentioned alternative approaches. 

As in refs. [10,11] we concentrate on continuous-time systems and start by 
observing that a DN-term modeled by (the increments of) a Wiener process 
Sw k (coupled through some constants Gjk)' 

dxj = Fj(x)dt + ^2 Gjk^Wk (1) 
k 

is basically different from white AN, superimposed to the time series {s n = 
s[x(nAt)]} after a clean (Gj k = 0) integration: 

S n = S n + 7] n . (2) 



This is seen in the typical case in which the "measurement function" s(x) is 
simply one of the coordinates, say Xj. In fact, while the white- noise process rj n 
of eq. (2) appears to be always discontinuous in time, the Sw k in eq. (1) lead 
to a continuous solution, as it can be seen from the increments in a time 5t: 

xj(t + 5t) = xj(t) + F,St + ]T G jklk {St)^ 2 + ■■■ (3) 

k 

The 7fc's are 5-correlated normal Gaussian random numbers, and in the limit 
St — > one has Xj(t + St) — > xj(t) [18]. The same applies for a generic (at 
least differentiable) measurement function, as can be seen from a first-order 
expansion: 

s[x(t + St)] = s[x(t)\ + £(^)[i^t + £ G jklk {Stf' 2 ] . (4) 

3 k 

Apart from this basic difference, the distinction between AN and DN in the 
real output of a numerical integration (or of an experimental device) is still a 
complicated task because the continuity must be judged on the basis of finite 
increments At over which the series acquires a finite increment As. From this 
point of view, the continuity statistics (CS) of Pecora et al. [11,13] sheds some 
light as it evaluates the degree of continuity at different resolution scales. 
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Before closing the section we draw some comments on another type of problem 
which is sometimes put forth [19,2], namely, that the distinction between AN 
and DN is not well-posed since (at least in some cases) they can be mapped 
onto each other. Consider a discrete-time dynamics with some noisy feedback 
f. 

x n+1 = f (x n ) + 7 n+ i , n = 0,l,... (5) 

In absence of DN the clean time series from the j-ih coordinate would be: 

SO = XOj , Si = fj (x ) , S 2 = fj (f(x )) , s 3 = fj (/ (/(f ))) , • • • (6) 

whereas DN turns it to: 

s = x 0j , si = fj (x ) + 7i,-, s 2 = fj (/ (f ) + 7i) + 72 j, 

s 3 = fj (/ (M) + 71) + 72) + 73, ... (7) 
Now, by the light of eq. (2), one could equally say that the time series (7) stems 

— # 

from a deterministic dynamics x n +i = f (x n ) plus the following (equivalent) 
AN: 

Vi = lij, V2 = l2j + fj (f(x ) + 7i) - fj (/ (^o)) , 
Vs = 73, + /, (/ (/ (fo) + 7i) + 72) - fi (f(f(xo))) (8) 

Similar arguments apply to more general measurement functions and also to 
continuous-time systems, at least at the leading order of eq. (4). Thus, in prin- 
ciple it is conceivable that a special type of AN can mimic a given DN present 
in the dynamics. However, from a practical point of view, a measurement noise 
like the one in eq. (8) is not too realistic, as it is "customised" on the dynamics 
itself (through /) and on the initial conditions (in the case of continuous-time 
systems it would also scale with the integration time). In addition, what is 
more important is that the process in eq. (8) is autocorrelated, as can be seen 
by inverting the relationship between the 7's and the 77's step by step. Again, 
the autocorrelation pattern is determined by dynamical features of the system 
in a complicated way and in general it does not correspond to a typical mea- 
surement component. Thus, being conscious of the mathematical subtleties 
which may arise in rigorous framework, in what follows we restrict to white 
noise processes as a first possible modelisation of real situations. 
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3 Model and Numerical Procedures 



3. 1 The Van der Pol oscillator 



The different effects of AN and DN have been tested on the Van der Pol 
system, whose "clean" equations are: 

y = -jj,( x 2 -l)y-x 

x = y (9) 



The white noise, added either to the left-hand sides of eqs. (9) or to the clean 
coordinate x, is generated using a Gaussian distribution with zero mean. Dif- 
ferent values of the variance, a 2 , have been considered in order to tune the 
strength of the noisy terms. We point out that from the discussion of the pre- 
vious section, and from the following results, there are no apparent obstacles 
to extend these ideas to high-dimensional/chaotic systems. The reason why 
we choose such a simple system is that, due to its one-dimensional attractor, 
we can clarify the points without worrying about the problems of "wandering" 
[11] which may arise when the embedding dimension is pushed at higher and 
higher values. 



3.2 Numerics 



The initial conditions for the system (9) have been set to (xq — l,yo — 0) 
and n — 1. These choices lead to a "clean" amplitude of oscillation of about 
2, with period 6.6. Dealing with stochastic differential equations, we used an 
Euler integration scheme to produce time-series of 16384 points from the x 
coordinate. To achieve a good compromise in terms of numerical accuracy we 
choosed an integration time step At = 0.01. The transient time was observed 
to be sufficient for the oscillator to relax on the ID attractor. Numerically we 
explored embedding dimension m up to 10, using a delay time of 166 which 
corresponds to the first zero of the autocorrelation function (in units of At). 
As in [11], we adopt an Epanechnikov kernel measure estimator with a radius 
of 0.05 (after that the reconstructed attractor has been rescaled to lie within 
the hypercube [0, 1] x • • • x [0, 1]). 
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3.3 Continuity statistics 



In order to test the mathematical properties embodied in a possible mapping 
between two given time series, that is, continuity, differentiability, inverse dif- 
ferentiability and injectivity, Pecora et al. [13] have developed a set of statistics 
aimed to test quantitatively these features. Their algorithms are of general use 
and can in particular be applied to test topological properties in any pair of 
sets of points. Basically, the method is intended to evaluate, in terms of prob- 
ability or confidence levels, whether two data sets are related by a mapping 
having the continuity property: A function / is said to be continuous at a 
point xo if Ve > 0, 35 > such that || x — x \\< 5 =>\\ f(x)) — f(x ) \\< e. The 
results are tested against the null-hypothesis, specifically, the case in which 
no functional relation exists. This is done by means of the statistics proposed 
by Pecora et al. [13] 

©co(e) = -Ee c a(e,j) (10) 
u p j=i 

and 



Qco(e,j) = l--^- (11) 

max 

where pj is the probability that all of n§ points in the 5- set, around a certain 
point Xj, fall in the e-set around f(xj). The likelihood that this will happen 
must be relative to the probability, P max , of the most likely event under the 
null hypothesis. In the Appendix we present some details on the calculation 
of the ratio in eq. (11), useful for a numerical implementation. The sum in eq. 
(10) represent an average over n p points chosen at random in the whole time 
series. Now, when B^o (e) » 1 we can confidently reject the null hypothesis, 
and assume that there exists a continuous function. As in the work of Pecora 
et al. [13] the e scale is relative to the standard deviation of the time series, 
and thus e G [0,1]. Plots of Q c ° (s) versus e can be used to quantify the 
degree of statistical continuity of a given function. In order to characterise the 
continuity statistics by means of a single parameter we have also calculated, 

i 

= J e c0 (e)de (12) 
o 



The limiting values of 6, namely, and 1, correspond to a strongly discon- 
tinuous and a fully continuous function, respectively. In the results reported 
here the function / mentioned above can be either the logarithmic derivative 
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of the measure or the coordinate itself. When useful, we will denote the cor- 
responding continuity statistics by CSLDM and CSC. The resolution scale e 
varied in the range [0.0001, 1] and n p was always fixed to be a 20% of the total 
length of the series. 



4 Results 

The CSC for the amplitude x of the Van der Pol oscillator (eqs. (9)) is shown 
in fig. 1 for different levels and kinds of noise. On the one hand it is readily seen 
that the CSC for the time series affected by DN is essentially independent of 
the noise amplitude. This appears to be consistent with the continuity analysis 
sketched in section 2. On the other hand, the CSC for series affected by AN 
decreases steadily, and almost linearly, with the noise level a an- We have 
considered also the more subtle case in which both AN and DN are present. 
In particular, we have fixed a moderate level of DN, like o"dn — 0.1, and 
then contaminated the resulting x-time series with various levels of AN. The 
corresponding points in fig. 1 are almost indistinguishable from those without 
DN, indicating that what rules the CSC is the presence of AN. To give an 
idea of how the CS is disributed over the different resolution scales we have 
plotted, in fig. 2, the Qc° statistics of the coordinate for some representative 
cases. Roughly we could say that the effect of the various noise sources is to 
shift, along the e axe, the same sigmoidal shape. What is different, between 
AN and DN, is that the shift induced by the former is considerably more 
pronounced, actually one order of magnitude in e when passing from a an = 

0.03 tO (T AN = 0.3. 

Let us now come to the analysis of the CSLDM values. In refs. [11,17] it was 
argued that this quantity is an indicator which is sensible to the presence of 
DN. The results of fig. 3 show that it is also strongly affected by AN. Both 
the DN- and the AN-data decrease almost monotonically with the noise level, 
with the exception of a small jump at a — 0.09. The origin of the latter is not 
completely clear, though it may be partially due to the statistical fluctuations 
induced by the sampling over the n p points (eq. (10)). As in fig. 1 we present 
some cases in which the combined action of AN and DN takes place, just to 
show a limitation of the present approach. The pure AN and the AN+DN data 
are not so close as in fig. 1 but there is not a definite trend to discriminate 
between a non-stochastic (odn = 0) and a stochastic underlying dynamics, 
when the original time series is contaminated by AN. As a general property 
we should note that the AN-free values are always greater than the ones in 
which AN is present. For the sake of completeness a O-vs-s plot for the measure 
is given in fig. 4. With respect to fig. 2 the sigmoidal shape is flattened (rather 
than shifted), and the AN curves are somehow less regular than the DN ones. 
Nonetheless, the lowering effect of noise is clearly appreciable. Note also that 
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in this case the effect of varying the AN level of one order of magnitude is far 
less evident than in the CSC. In fact this is consistent with fig. 3 where it is 
seen that the CSLDM jumps down and levels off to ~ 0.3 as soon as a small 
amount of AN is introduced. 



5 Conclusions 

In this paper we have presented the extension of the method of refs. [10,11] 
to the analysis of time series affected by AN. The basic task which one has 
to perform is to compare the behaviour of the statistical continuity of the 
coordinate and of the statistical differentiabilty of the natural measure, at 
different levels of noise. While the first is sensibly different from the "clean" 
one only when AN is present, the second is affected by both dynamical and 
additive noise. Hence, it is possible to discriminate between the two cases. 
Altough we believe that the present method can be readily applied also to 
high-dimensional and/or chaotic systems, the results discussed here refer to the 
simple Van der Pol oscillator. A further key step in the analysis of experimental 
time series, namely the criterion to adopt when both types of noise are present, 
has been only touched here and is currently investigated by means of real 
physiological data. 
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7 Appendix 

As pointed out in [13] the appropriate probability distribution for the conti- 
nuity test is the binomial one: 

P(k; P, n s ) = kl (" s S l k y P k (l - P) ns ~ k , k = 0, 1, . . . , n s , (13) 

with p = n £ /N, n £ being the number of points in e-set (see text). In addition, 
due to the null hypotesis under consideration, the probility pj appearing in 
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eq. (11) is just pj = p n& . Except for the trivial case p — 1, for which one 
has P(k) = 5k,n S i the maximum of the distribution (13) is located at /c max — 
[p(n<5 + 1)] ([ • ] denoting the integer part - see, for istance, [20]). Now, let us 
observe that in the present case we must calculate only the ratio: 

P H6 k max \ (ns — fc max) ■ -~m max 



P max ^<5 ■ 



where p = p/(l — p). /,From a numerical point of view it is advantageous 
to exploit the following trick in the left-hand side of eq. (14). Take M = 
max (k max , ns — A; max ), so that Ml can be simplified in the numerator and in 
the denominator of the fraction leaving: 

k m J-(n 5 - fc max )! _ (n s - M)\ _ i 

n & \ n s (n 5 - 1) . . . (M + 1) ^ t + M ' 1 j 



Dealing with the product in eq. (15) has the advantage of avoiding ratios 
of very large numbers. However, the powers of p appearing in eq. (14) can 
still introduce rather small numbers, which have to be treated through their 
logarithms. 
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Fig. 1. CSC versus noise standard deviation, a, for the x coordinate of the Van der 
Pol oscillator. Results are shown for pure AN (no DN, open circles), pure DN (no 
AN, full diamonds) and for a system having a fixed level ctdn = 0-1 plus AN with 
variable standard deviation (full squares, almost overlapped with the open circles). 




Fig. 2. CSC resolved at different e scales. The legend indicates the type and the 
levels (standard deviation) of noise. 
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Fig. 3. CSLDM versus noise standard deviation, a, corresponding to a 
four-dimensional embedding (the other paramters being reported in the main text). 
Results are shown for pure AN (no DN, open circles), pure DN (no AN, full dia- 
monds) and for a system having a fixed level <tdn = 0.1 plus an AN with variable 
standard deviation (full squares). 
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Fig. 4. CSLDM resolved at different e scales, for the same levels and types of noise 
of Fig. 2. In all cases, the measure is estimated from a four-dimensional embedding 
with the parameters indicated in the text. 
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